package jpave.qr;

import java.awt.FileDialog;
import java.awt.Frame;
import java.io.BufferedReader;
import java.io.InputStreamReader;
import java.sql.Connection;
import java.sql.PreparedStatement;
import java.sql.ResultSet;
import java.sql.Statement;
import java.util.Collection;
import java.util.Iterator;
import java.util.TreeMap;
import java.util.TreeSet;
import java.util.Vector;
import java.util.HashSet;
import java.util.Map;
import java.lang.StringBuilder;

import org.rosuda.JRI.REXP;
import org.rosuda.JRI.RFactor;
import org.rosuda.JRI.RList;
import org.rosuda.JRI.RMainLoopCallbacks;
import org.rosuda.JRI.RVector;
import org.rosuda.JRI.Rengine;

import assem.Schema;

import util.Debug.ErrorReport;
import util.methods.TimeHelpers;
import jpave.dataholders.*;
import assem.DBConn;

public class QRProcess {
	public static final int METHOD_EDGER = 1;
	public static final int METHOD_DEGSEQ = 2;
	public static final int METHOD_BAYSEQ = 3;
	public static final int METHOD_DESEQ = 4;
	
	public static final int NORM_DEFAULT = 1;
	public static final int NORM_EDASEQ = 2;
	public static final int NORM_NONE = 3;
	
	public static final int MINOR_METHOD_ORIGINAL = 1;
	public static final int MINOR_METHOD_GLM = 2;
	public static final int MINOR_METHOD_LRT = 3;
	public static final int MINOR_METHOD_CTR = 4;
	public static final int MINOR_METHOD_FET = 5;
	public static final int MINOR_METHOD_MARS = 6;
	public static final int MINOR_METHOD_FC = 7;
	
	public static final String [] METHOD_NAMES = { "edgeR", "DEGseq", "DESEq"};
	public static final int [] METHOD_IDS = { METHOD_EDGER, METHOD_DEGSEQ, METHOD_DESEQ};
	public static final int [] NORM_IDS = { NORM_DEFAULT, NORM_EDASEQ };
	public static final String [] NORM_NAMES = { "Default", "EDAseq" };
	
	public static final int [] EDGER_NORM_IDS = { NORM_DEFAULT, NORM_EDASEQ };	
	public static final String [] EDGER_MINOR_NAMES = {};// "original", "GLM" };
	public static final int [] EDGER_MINOR_IDS = {}; // MINOR_METHOD_ORIGINAL, MINOR_METHOD_GLM };

	public static final String [] DESEQ_MINOR_NAMES = {};
	public static final int [] DESEQ_MINOR_IDS = {}; 

	public static final String [] DEGSEQ_MINOR_NAMES = { "LRT", "FET", "MARS" };
	public static final int [] DEGSEQ_MINOR_IDS = { MINOR_METHOD_LRT, MINOR_METHOD_FET, MINOR_METHOD_MARS};
	
	public static final String [] BAYSEQ_MINOR_NAMES = {};
	public static final int [] BAYSEQ_MINOR_IDS = {};
	
	public static double dispersion = 0.1;
	
	private static String doPlots = "";

	public QRProcess(SessionData s, Connection c) {
		selectedDB = s;
		conn = c;
	}

	public void doExecute(int method, int normalization, int MINOR_METHOD,
		String pColName, boolean addCol, boolean doFDR, TreeSet<String> g1, TreeSet <String> g2, double disp) 
	{
		try
		{
			Statement s = selectedDB.getConnection().createStatement();
			ResultSet rs;
			String sql = "select schemver from schemver where schemver='" + Schema.currentVerString() + "'";
			rs = s.executeQuery(sql);
			if (!rs.first())
			{
				System.err.println("Database is not current, view database with viewSingleTCW to update");
				System.err.println("Terminate executing DE");
				return;				
			}
			TreeSet<String> grp1 = g1;
			TreeSet<String> grp2 = g2;	

			long startTime = System.currentTimeMillis();
			
			initJRI();
			packages = new HashSet<String>();
			getPackages(packages);
			
			int[] samps = countsToR(grp1,grp2);
			int nSamps1 = samps[0];
			int nSamps2 = samps[1];
			
			if (normalization == NORM_EDASEQ)
			{
				if (!checkPackage("EDASeq")) return;
				EDANorm();
			}
			TreeMap<String,Double> scores = new TreeMap<String,Double>();
			if (method == METHOD_EDGER)
			{
				runEdgeR(scores, normalization,nSamps1,nSamps2, disp);
			}
			else if (method == METHOD_DEGSEQ)
			{
				runDEGseq(scores,normalization,MINOR_METHOD);
			}
			else if (method == METHOD_DESEQ)
			{
				runDESeq(scores,normalization,nSamps1,nSamps2, disp);
			}				
			else
			{
				System.err.println("Error in method selected");
				return;
			}

			if (scores.size() > 0 && doFDR)
			{
				doFDR(scores);
			}
			if (scores.size() > 0 && pColName.length() > 0 && addCol)
			{
				addDEcols(scores, pColName, grp1, grp2);
			}       
			System.err.println("\nFinished DE execution at: " + TimeHelpers.getElapsedTimeStr(startTime));
			System.err.println("\n" + doPlots);
        		re.startMainLoop();
			System.err.println("The console is in R, you may run R commands -- quit() or Cntl-C " +
					"when done, or perform another Execute.");
			System.err.print(">");
		}
		catch(Exception e)
		{
			ErrorReport.reportError(e, "Error in execute");
		}
	}
	private void addDEcols(TreeMap<String,Double> scores, String pColName, 
			TreeSet <String> grp1, TreeSet <String> grp2) 
	{
		try {
			System.err.println("\nLoading " + scores.size() + " scores");
			Statement s = selectedDB.getConnection().createStatement();
			ResultSet rs = s.executeQuery("show columns from contig where field= '" + pColName + "'");
			if (!rs.first())
			{
				String tmp = pColName.substring(2); // so doesn't print P_
				System.err.println("   Adding database column " + tmp + " (may take a while)...");
				s.executeUpdate("alter table contig add " + pColName + " double default 2");
			}
			else
			{
				s.executeUpdate("update contig set " + pColName + "=2");
			}
			s.executeUpdate("update assem_msg set pja_msg=NULL");
		
			PreparedStatement ps = conn.prepareStatement("update contig set " + pColName + "=? where contigid=?");
			int num = 0;
			for (String ctg : scores.keySet())
			{
				if (scores.get(ctg).isNaN()) continue;
				ps.setDouble(1,scores.get(ctg).doubleValue());
				ps.setString(2, ctg);
				ps.addBatch();
				num++;
				if (num % 100 == 0) 
				{
					System.err.print(num + "      \r");
					try
					{
						ps.executeBatch();
					}
					catch(Exception e)
					{
						System.err.println("BAD VALUE " + ctg + " " + scores.get(ctg));
					}
				}
			}
			if (num % 100 != 0) ps.executeBatch();
			
			// Now IF there is one library in each group set, and IF the corresponding ln__ columns exist, 
			// then we will convert the scores to +/- to encode the direction of DE.
			if (grp1.size() == 1 && grp2.size() == 1)
			{
				String col1 = "LN__" + grp1.first();
				String col2 = "LN__" + grp2.first();
				rs = s.executeQuery("show columns from contig where field='" + col1 + "'");
				boolean col1Exists = rs.first();
				rs = s.executeQuery("show columns from contig where field='" + col2 + "'");
				boolean col2Exists = rs.first();
				if (col1Exists && col2Exists)
				{
					s.executeUpdate("update contig set " + pColName + " =-" + pColName + 
							" where " + col1 + "<" + col2);						
				}
			}
			// XXX CAS 12/28/14 add LibraryDE.title
			String title = "";
			for (String lib : grp1) title += lib + " ";
			title += ": ";
			for (String lib : grp2) title += lib + " ";
			if (title.length()>100) title = title.substring(0,99);
		
			rs = s.executeQuery("show tables like 'libraryDE'");
			if (!rs.first()) {
				System.out.println("   Create library DE table");
				s.executeUpdate("create table libraryDE (" +
						"pCol varchar(30), title varchar(100), " +
						"index(pCol));");
			}
			rs = s.executeQuery("Select title from libraryDE where pCol='" + pColName + "'");
			if (rs.first()) {
				s.executeUpdate("update libraryDE set title='" + title + 
						"' where pCol='" + pColName + "'");
			}
			else {
				s.executeUpdate("insert into libraryDE set pCol='" + 
						pColName + "', title='" + title + "'");
			}
			System.err.println("   Finish adding scores to sTCW");	
		}
		catch (Exception e) {ErrorReport.prtReport(e, "Adding DE columns");}
	}
	void doFDR(TreeMap<String,Double> ctg2score )
	{
		Vector<String> names = new Vector<String>();
		double[] scores = new double[ctg2score.size()];
		int i = 0;
		for (String name : ctg2score.keySet())
		{
			names.add(name);
			scores[i] = ctg2score.get(name);
			i++;
		}
		try
		{
			if (re == null) initJRI();
        	String cmd;
        	REXP x;
        	
        	System.err.println("\nConvert to FDR....");
        	
        	cmd = "library(\"multtest\")";
        	System.err.println("     " + cmd);
        	re.eval(cmd);

        	re.assign("rawScores", scores);
			
        	cmd = "scoresFDR <- mt.rawp2adjp(rawScores,c(\"BH\"))";
        	System.err.println("     " + cmd);
        	re.eval(cmd);
        	
        	cmd = "orderedFDR <- scoresFDR$adjp[order(scoresFDR$index),2]";
        	System.err.println("     " + cmd);
        	x = re.eval(cmd);
        	
        	//cmd = "finalFDR <- orderedFDR[,2]";
        	//System.err.println(cmd);
        	//x = re.eval(cmd);        	
        	
        	double[] adjScores = x.asDoubleArray();
        	assert(adjScores.length == names.size());
        	for (i = 0; i < names.size(); i++)
        	{
        		String name = names.get(i);
        		double fdr = adjScores[i];
        		ctg2score.put(name, fdr);
        	}
        	System.err.println("FDR done");
		}
		catch(Exception e)
		{
			ErrorReport.reportError(e, "Error computing FDR");
		}
	}
	void getPackages(HashSet<String> pkgs)
	{
	    	REXP x = re.eval("rownames(installed.packages())");
	    	for (String pkg : x.asStringArray())
	    	{
	    		pkgs.add(pkg);
	    	}		
	}
	boolean checkPackage(String pkg)
	{
		if (!packages.contains(pkg)) 
		{
			System.err.println("*****R package " + pkg + " needs to be installed.******");
			return false;
		}
		return true;
	}
	/***
	 *  Put the expression data into R objects for use by the DE methods
	 */
	int[] countsToR(TreeSet<String> grp1, TreeSet<String> grp2)
	{
		try
		{
			Statement s = selectedDB.getConnection().createStatement();
			ResultSet rs;
/* 			backwards compatible for databases that do not have the rep field 
*			Now should be caught by DB check, unless version is somehow 4.0 but fields still missing. 				
*/
			rs = s.executeQuery("show columns from library where field='reps'");
			if (!rs.first()) {
        			s.executeUpdate("alter table library add reps text after parent");
        			s.executeUpdate("update library set reps=''"); // not tested, did intereactively
			}
			rs = s.executeQuery("show columns from clone_exp where field='rep'");
			if (!rs.first())
			{
				System.err.println("Add rep column to clone_exp...");
				s.executeUpdate("alter table clone_exp add rep smallint unsigned default 0");
			}
			
			Vector<String> allGrp = new Vector<String>();
			allGrp.addAll(grp1);allGrp.addAll(grp2);
			Vector<String> allGrpReps = new Vector<String>(); // grp names, or expanded to rep names		
			
			// Assign each lib or replicate an order
			TreeMap<Integer,TreeMap<Integer,Integer>> sampOrder = new TreeMap<Integer,TreeMap<Integer,Integer>>();
			int totalSamps = 0;
			int nSamps1 = 0; int nSamps2 = 0;
			s.executeUpdate("update library set reps='' where reps is null"); // in case...
			System.err.println("\n"); // makes replicate lines slightly prettier, if any
			for (String lib : allGrp)
			{
				rs = s.executeQuery("select reps,lid from library where libid='" + lib + "'");
				rs.first();
				String reps = rs.getString(1).trim();
				int lid = rs.getInt(2);
				sampOrder.put(lid,new TreeMap<Integer,Integer>());
				if (reps.equals(""))
				{
					sampOrder.get(lid).put(0, totalSamps);
					totalSamps++;
					if (grp1.contains(lib)) nSamps1++;
					else nSamps2++;
					allGrpReps.add(lib);
				}
				else
				{
					String[] repnames = reps.split(",");
					int nReps = repnames.length;
					System.err.println(lib + " has " + nReps + " replicates:" + reps);
					for (int r = 1; r <= repnames.length; r++)
					{
						sampOrder.get(lid).put(r, totalSamps);
						totalSamps++;	
						if (grp1.contains(lib)) nSamps1++;
						else nSamps2++;
						allGrpReps.add(repnames[r-1].trim());
					}
				}
			}
			System.err.println("Collecting count data (may take several minutes)");
			
			rs = s.executeQuery("select count(*) from contig");
			rs.first();
			int nCtg = rs.getInt(1);
			String[] names = new String[nCtg];
			double[] gc = new double[nCtg];
			
			// Don't change the contig order in next query as GOseq depends on it
			rs = s.executeQuery("select contigid, ctgid, gc_ratio from contig order by ctgid asc limit " + nCtg); // for testing
			int ctgNum = 0;
			TreeMap<Integer,Integer> ctg2num = new TreeMap<Integer,Integer>();
			while (rs.next())
			{
				names[ctgNum] = rs.getString(1);
				gc[ctgNum] = rs.getFloat(3);
				ctg2num.put(rs.getInt(2), ctgNum);
				ctgNum++;
			}
			int[] cnts = new int[totalSamps*nCtg];
			
			String ctgClause = "";
			String libClause = " and lid in (" + strJoin(sampOrder.keySet(),",") + ") ";
			rs = s.executeQuery("select ctgid, lid, rep,sum(count) from clone_exp join contclone " +
						" on contclone.cid=clone_exp.cid where 1 " + ctgClause + libClause + " group by contclone.ctgid, lid, rep");							
			while (rs.next())
			{
				int ctgid = rs.getInt(1);
				int cnum = ctg2num.get(ctgid);
				int lid = rs.getInt(2);
				int rep = rs.getInt(3);
				int snum = 0;
				if (sampOrder.get(lid).containsKey(0))
				{
					snum = sampOrder.get(lid).get(0);
				}
				else if (sampOrder.get(lid).containsKey(rep))
				{
					snum = sampOrder.get(lid).get(rep);
				}
				else
				{
					continue;
					// It could be the total count (rep=0) for a library that has reps
				}
				cnts[cnum + snum*nCtg] = rs.getInt(4);
			}
			if (re == null) initJRI();
     
        	String cmd;
        	
        	assert(allGrpReps.size() == nSamps1 + nSamps2);
        	System.err.println("Total samples " + totalSamps);
        	re.assign("numbers", new int[]{nSamps1,nSamps2,totalSamps,nCtg});
        	cmd = "nSamps1 <- " + nSamps1;
        	System.err.println("     " + cmd);
        	re.eval(cmd);
        	cmd = "nSamps2 <- " + nSamps2;
        	System.err.println("     " + cmd);
        	re.eval(cmd);
        	cmd = "totalSamps <- " + totalSamps;
        	System.err.println("     " + cmd);
        	re.eval(cmd);
        	cmd = "nCtg <- " + nCtg;
        	System.err.println("     " + cmd);
        	re.eval(cmd);

        	re.assign("names",names);
        	re.assign("gc",gc);
        	re.assign("grpnames",allGrpReps.toArray(new String[0]));
        	re.assign("cnts", cnts);
        	cmd = "c <- array(cnts,dim=c(" + nCtg + "," + totalSamps + "))";
        	System.err.println("     " + cmd);
        	re.eval(cmd);
        	cmd = "rownames(c) <- names";
        	System.err.println("     " + cmd);
        	re.eval(cmd);  
        	
        	return new int[]{nSamps1,nSamps2};
		}
		catch(Exception e)
		{
			ErrorReport.die(e, "performing DB counts");
		}
		return null;
	}
	
	private String strJoin(Collection s, String delimiter)
	{
		StringBuilder buffer = new StringBuilder();
		Iterator iter = s.iterator();
		while (iter.hasNext())
		{
			buffer.append(iter.next().toString());
			if (iter.hasNext())
				buffer.append(delimiter);
		}
		return buffer.toString();
	}
	void finishGC()
	{
		try
		{
			Statement s = selectedDB.getConnection().createStatement();
			ResultSet rs;
		
			rs = s.executeQuery("select count(*) from contig where gc_ratio is null");
			rs.first();
			if (rs.getInt(1) > 0)
			{
				System.err.println("Complete gc ratio calculations");
				s.executeUpdate(" update contig set gc_ratio = " + 
						" length(replace(replace(replace(replace(replace(replace(consensus,\"n\",\"\"),\"N\",\"\"),\"t\",\"\"),\"a\",\"\"),\"T\",\"\"),\"A\",\"\"))/" +
						"length(replace(replace(consensus,\"n\",\"\"),\"N\",\"\")) " +
						"  where gc_ratio is null and consensus_bases > 0");
						
				s.executeUpdate("update contig set gc_ratio=0 where gc_ratio is null");
			}
		}
		catch(Exception e)
		{
			ErrorReport.reportError(e, "Computing GC");
			return;
		}
		
	}
	/**********************************************************
	 * EDANorm
	 */ 
	void EDANorm()
	{
		if (!checkPackage("EDASeq")) return;

		finishGC(); 

		System.err.println("\nBegin EDA normalization");
		re.eval("library(EDASeq)");
		
		String cmd;
	    	cmd = "phenoData <- data.frame(conditions=c(rep(\"grp1\",nSamps1),rep(\"grp2\",nSamps2)),row.names=grpnames)";
	    	System.err.println("     " + cmd);
	    	re.eval(cmd);

	    	cmd = "feature <- data.frame(gc=gc,row.names=names)";
	    	System.err.println("     " + cmd);
	    	re.eval(cmd);

	    	cmd = "myData <- newSeqExpressionSet(exprs=c, featureData=feature, phenoData=phenoData)";
      	System.err.println("     " + cmd);
      	re.eval(cmd);
    	
      	cmd = "dataWithin <- withinLaneNormalization(myData,\"gc\",which=\"full\")";
      	System.err.println("     " + cmd);
      	re.eval(cmd);

      	cmd = "dataNorm <- betweenLaneNormalization(dataWithin,which=\"full\")";
      	System.err.println("     " + cmd);
      	re.eval(cmd);  	
	}
	/*******************************************************
	 * runEdgeR
	 */
	private void runDESeq(TreeMap<String,Double> scores,int normalization,int nSamps1, int nSamps2, double disp)
	{
		if (!checkPackage("DESeq")) return;
		System.err.println("\nBegin DESeq");
		
		try
		{		
        	REXP x;
        	String cmd;
        	re.eval("library(DESeq)");
        	
        	cmd = "group <- c(rep(\"Grp1\",nSamps1),rep(\"Grp2\",nSamps2))"; 
        	System.err.println("     " + cmd);        	
        	re.eval(cmd);
        	      	
        	if (normalization == NORM_EDASEQ)
        	{
            	cmd = "cds <- newCountDataSet(exprs(dataNorm),group)";
        		System.err.println("     " + cmd);
        		re.eval(cmd);
        	}
        	else
        	{
            	cmd = "cds <- newCountDataSet(c,group)";
        		System.err.println("     " + cmd);
        		re.eval(cmd);
        	}
        	
        	cmd = "cds <- estimateSizeFactors(cds)";
    		System.err.println("     " + cmd);
    		re.eval(cmd);

    		// The next command has to be run before manually setting dispersions, to create the structure
			// fitType is needed else it sometimes doesn't converge
    		String estCmd = "cds <- estimateDispersions(cds,method='blind',sharingMode='fit-only',fitType='local')";

    		if (disp >= 0) // this means it was set by user
    		{
	        	cmd = estCmd;
	    		System.err.println("     " + cmd);
	    		re.eval(cmd);
	    		
	        	cmd = "fData(cds)[,1] = rep(" + disp + ",length(fData(cds)[,1]))";
	    		System.err.println("     " + cmd);
	    		re.eval(cmd);
    		}
    		else if (nSamps1 > 1 || nSamps2 > 1)
    		{
	        	cmd = "cds <- estimateDispersions(cds,sharingMode='gene-est-only')";
	    		System.err.println("     " + cmd);
	    		re.eval(cmd);
    		}
    		else
    		{
    			System.err.println("*** No replicates! Dispersion estimate may be inaccurate. ***");
    			System.err.println("*** It is recommended to set dispersion manually. ***");
	        	cmd = estCmd;
	    		System.err.println("     " + cmd);
	    		re.eval(cmd);	    		
    		}

        	cmd = "res <- nbinomTest(cds,\"Grp1\",\"Grp2\")";
    		System.err.println("     " + cmd);
    		re.eval(cmd);
    		
        	cmd = "res$pval[is.na(res$pval)] <- 2";
    		System.err.println("     " + cmd);
    		re.eval(cmd);
  		
    		RList results = re.eval("res").asList();
       	 
        	String[] ctgs = results.at("id").asStringArray();
        	double[] pvals = results.at("pval").asDoubleArray();

        	for (int i = 0; i < ctgs.length; i++)
        	{
        		scores.put(ctgs[i], pvals[i]);
        	}

        	System.err.println("DESeq done");
        	doPlots = "DESeq plots:\n" +
        			  "  Plot dispersion esimates: plotDispEsts( cds )\n" +
        			  "  MA plot: plotMA(res)";     
        }
        catch(Exception e)
        {
	        	ErrorReport.die(e, "Executing DESeq");
        }		
	}
	
	/*******************************************************
	 * runEdgeR
	 */
	private void runEdgeR(TreeMap<String,Double> scores,int normalization,int nSamps1, int nSamps2, double disp)
	{
		if (!checkPackage("edgeR")) {
			System.err.println("EdgeR package does not exist");
			return;
		}
		System.err.println("\nBegin EdgeR");
		
		try
		{		
	        	REXP x;
	        	String cmd;
	        	re.eval("library(edgeR)");
	        	
	        	String grpStr = "group <- c(rep(\"Grp1\",nSamps1),rep(\"Grp2\",nSamps2))"; 
	        	re.eval(grpStr);
	        	      	
	        	if (normalization == NORM_EDASEQ)
	        	{
	        		cmd = "d <- DGEList(counts = exprs(dataNorm), group = group)";
	        		System.err.println("     " + cmd);
	        		re.eval(cmd);
	        	}
	        	else
	        	{
	        		cmd = "d <- DGEList(counts = c, group = group)";
	        		re.eval(cmd);
	                // Normalization: p10 using the default normalizing method TMM 	
	                // computes effective library size -  a set of scaling factors for the library sizes 
	               	// that minimize the log-fold changes between the samples for most genes.
	
	        		cmd = "d <- calcNormFactors(d)";
	            	System.err.println("     " + cmd + " (using edgeR default TMM normalization)");
	            	re.eval(cmd);
	        	}
	        	
	        	if (disp >= 0) // this means it was set by user
	        	{
	        		cmd = "d$common.dispersion <- " + disp;
		        	System.err.println("     " + cmd);
		        	re.eval(cmd);	        		
	        	}
	        	else if ((nSamps1 > 1 || nSamps2 > 1))
	    		{
		        	cmd = "d <- estimateCommonDisp(d)";
		        	System.err.println("     " + cmd);
		        	re.eval(cmd);

		        	cmd = "d <- estimateTagwiseDisp(d)";
		        	System.err.println("     " + cmd);
		        	re.eval(cmd);
	    		}
	    		else
	    		{
	    			disp = QRProcess.dispersion;
	    			System.err.println("*** No replicates! Assuming dispersion=" + disp + " ***");
	        		cmd = "d$common.dispersion <- " + disp;
		        	System.err.println("     " + cmd);
		        	re.eval(cmd);
	    		}

	        	// tests for DE
	        	cmd = "et <- exactTest(d)$table";
	        	System.err.println("     " + cmd);
	        	re.eval(cmd);
	        	
	        	cmd = "et$PValue[is.na(et$PValue)] <- 2";
	        	System.err.println("     " + cmd);
	        	re.eval(cmd);
	        	
	        	if (re.eval("et")==null) {
	        		System.err.println("Error: Cannot compute edgeR for these libraries");
	        		return;
	        	}
	        	RList results = re.eval("et").asList();
	 
	        	cmd = "rownames(et)";
	        	System.err.println("     " + cmd);
	        	x = re.eval(cmd);
	        	
	        	String[] ctgs = x.asStringArray();
	        	double[] pvals = results.at("PValue").asDoubleArray();
	
	        	for (int i = 0; i < ctgs.length; i++)
	        	{
	        		scores.put(ctgs[i], pvals[i]);
	        	}
	
	        	System.err.println("EdgeR done");
	        	doPlots = "EdgeR plots:\n" +
	        			  "  Scatter plot: plotBCV(d)\n" + 
	        			  "  Cluster plot: plotMDS(d, labels=grpnames)";     
        	}
        catch(Exception e)
        {
          	ErrorReport.die(e, "Executing edgeR");
        }	
	}
	/*********************************************************
	 * runDEGseq
	 */
	private void runDEGseq(TreeMap<String,Double> scores, int normalization, int methnum)
	{
		if (!checkPackage("DEGseq")) {
			System.err.println("DEGseq package does not exist");
			return;
		}
		System.err.println("\nBegin DEGseq");

    		String method = "";
		switch(methnum)
		{
			case MINOR_METHOD_LRT:
				method = "LRT"; break;
			//case MINOR_METHOD_CTR:
				//method =  "CTR"; break;
			case MINOR_METHOD_FET:
				method =  "FET"; break;
			case MINOR_METHOD_MARS:
				method =  "MARS"; break;
			//case MINOR_METHOD_FC:
				//method =  "FC"; break;
			default:
				System.err.println("DEGSeq bad method number " + method);
				return;
		}

		try
		{
	        	re.eval("library(DEGseq)");	
	        	REXP x;
	        	String cmd;
	        	        	
	        	//cmd = "outputDir <- file.path(\"DEGseqDir\", \"PAVE\")";
	        	cmd = "outputDir <- file.path(tempdir())";
	        	System.err.println("     " + cmd);
	        	x = re.eval(cmd);
	        	System.err.println("        type:" + x.getType());
	        	
	        	cmd = "m <- array(c(names,cnts),c(nCtg,totalSamps+1))";
	        	System.err.println("     " + cmd);
	        	re.eval(cmd);
	        	
	        	cmd = "DEGexp(geneExpMatrix1=m, geneCol1=1, expCol1=c(2:1+nSamps1), groupLabel1=\"grp1\", " +
	        			"geneExpMatrix2=m, geneCol2=1, expCol2=c((2+nSamps1):(totalSamps+1)), groupLabel2=\"grp2\","+
	        			"method=\"" + method + "\", outputDir=outputDir)";
	        	System.err.println("     " + cmd);
	        	re.eval(cmd);
	        	
	        	cmd = "resultFile <- file.path(outputDir, \"output_score.txt\")";
	        	System.err.println("     " + cmd);
	        	x = re.eval(cmd);
	
	        	cmd = "results <- read.table(resultFile,header=TRUE)";
	        	System.err.println("     " + cmd);
	        	re.eval(cmd);
	        	
	        	cmd = "results$q.value.Benjamini.et.al..1995.[is.na(results$q.value.Benjamini.et.al..1995.)] <- 2";
	        	System.err.println("     " + cmd);
        	
	        	x = re.eval("results");
	        	
	        	RVector results = re.eval("results").asVector();
	        	 
	        	RFactor genes = results.at("GeneNames").asFactor();
	
	        	double[] pvals = results.at("q.value.Benjamini.et.al..1995.").asDoubleArray();
	
	        	for (int i = 0; i < genes.size(); i++)
	        	{
	        		scores.put(genes.at(i), pvals[i]);
	        	}  	
     	   	System.err.println("DEGseq done");     
		}
		catch(Exception e)
		{
			ErrorReport.reportError(e, "running DEGseq");
		}
	}

	/**********************************************************
	 * runBaySeq (plus EDA)
	 */
/*	private void runBaySeq(TreeMap<String,Double> scores, int norm, int method)
	{
		if (!checkPackage("baySeq")) return;
		if (!checkPackage("snow")) return;

		try
		{
	        	REXP x;
	        	String cmd;
	        	re.eval("library(baySeq)");
	        	re.eval("library(snow)");
	        	
	    		cmd = "cl <- makeCluster(" + numThreads + ",\"SOCK\")";
	    		System.err.println("     " + cmd);
	    		re.eval(cmd);    		
	
	    		cmd = "replicates <- c(rep(1,nSamps1),rep(2,nSamps2))";
	    		System.err.println("     " + cmd);
	    		re.eval(cmd);
	
	    		// NDE model -- all groups
	        	// DE model -- all reps in first group is 1, all reps in second group is 2
	        	// if multiple libraries in either group, all reps are treated equally in the group
	        	cmd = "groups <- list(NDE = c(rep(1,totalSamps)) ,DE = c(rep(1,nSamps1),rep(2,nSamps2)))";
	        	System.err.println("     " + cmd);
	        	re.eval(cmd);
	        	
	        	if (norm == NORM_EDASEQ)
	        	{
	        		cmd = "counts <- exprs(dataNorm)";
	    			System.err.println("     " + cmd);
		        	re.eval(cmd);
		        	
				cmd = "libsizes <- colSums(counts)";
				System.err.println("     " + cmd);
				re.eval(cmd);
		
				cmd = "CD <- new(\"countData\", data=counts, replicates=replicates,libsizes=libsizes,groups=groups)";
				System.err.println("     " + cmd);
				re.eval(cmd);
	        	}
	        	else
	        	{
				cmd = "libsizes <- colSums(c)";
				System.err.println("     " + cmd);
				re.eval(cmd);
		
				cmd = "CD <- new(\"countData\", data=c, replicates=replicates,libsizes=libsizes,groups=groups)";
				System.err.println("     " + cmd);
				re.eval(cmd);  		
	        	}
			
			// baySEQ: negative binomial approach
			cmd = "CD <- getPriors.NB(CD,samplesize=10000, estimation=\"QL\",cl=cl)";
			System.err.println("     " + cmd);
			System.err.println("This takes awhile, even for a small dataset....");
			re.eval(cmd);
	
			cmd = "CD <- getLikelihoods.NB(CD, pET = \"BIC\", cl = cl)";
			System.err.println("     " + cmd);
			re.eval(cmd);
       	
	        	cmd = "tops <- topCounts(CD, group = \"DE\", number = nCtg)";
	        	System.err.println("     " + cmd);
	        	RList results = re.eval(cmd).asList();
	        	double[] pvals = results.at("FDR").asDoubleArray();
	        	 
	        	cmd = "rownames(tops)";
	        	System.err.println("     " + cmd);
	        	x = re.eval(cmd);
	        	String[] ctgs = x.asStringArray();
	
	        	for (int i = 0; i < ctgs.length; i++)
	        	{
	        		scores.put(ctgs[i], pvals[i]);
	        	}
	        	System.err.println("baySeq done");
	        	System.err.println("  MA plot: " +
	        			"plotMA.CD(CD, samplesA = \"1\", samplesB = \"2\")");
	     	System.err.println("  Posterior plot: " +
	     	  	"plotPosteriors(CD, group = \"DE\");");
		}
		catch(Exception e)
		{
			ErrorReport.reportError(e, "running baySeq -- try again.");
			return;
		}
	}
*/	
	/***********************************************
	 * GOSeq
	 */
	void runGOSeq(String[] colNames, String col, boolean doAll, boolean usePercent, double thresh)
	{
		String[] cols2do = (doAll ? colNames : new String[]{col});
		long startTime = System.currentTimeMillis();
		try
		{	    		
	        	ResultSet rs;
	        	String cmd;
	        int n = 0;
	       
	        	Statement s = selectedDB.getConnection().createStatement();
	        	rs = s.executeQuery("show tables like 'pja_unitrans_go'");
	    		if (!rs.first()) { // CAS 1/25/13
	    			System.err.println("\nWarning: The GO tables have not been added");
	    			System.err.println("Abort GOseq");
	    			s.close(); rs.close();
	    			return;
	    		}
	    		
	    		REXP x;
	        	initJRI();
	        	packages = new HashSet<String>(); // CAS added 1/25/13
	    		getPackages(packages);
	    		if (!checkPackage("goseq")) {
	    			System.err.println("GOseq package does not exist");
	    			return;
	    		}
	    		System.err.println("\nBegin GOseq");
	    	 	re.eval("library(goseq)");
	    		
	        	rs = s.executeQuery("select count(*) from contig where pid > 0");
	        	rs.first();
	        	int nCtg = rs.getInt(1);

	        	int[] lens = new int[nCtg]; 
	        	String[] names = new String[nCtg]; 
	        	rs = s.executeQuery("select contigid, consensus_bases from contig where pid > 0 order by ctgid asc");
	        	n = 0;
	        	while (rs.next())
	        	{
	        		String name = rs.getString(1);
	        		names[n] = name;
	        		int len = rs.getInt(2);
	        		lens[n] = len;
	        		n++;
	        	}
	        	re.assign("names",names);
	        	re.assign("ctgLens",lens);
			cmd = "nCtg <- " + nCtg;
			System.err.println("     " + cmd);
			re.eval(cmd);
				
	        	for (String colName : cols2do)
	        	{
	        		if (colName.equals("")) continue; // CAS 1/25/13
	        		if (colName.equals(QRFrame.selCol) || colName.equals(QRFrame.allCols)) continue;
	        		
	        		System.err.println("GOSeq: processing " + colName.substring(2));
	        		
				if (usePercent)
				{
					// get the percentage threshold
					if (thresh >= 100 || thresh <0)
					{
						System.err.println("Invalid percentage threshold " + thresh);
						return;
					}
					int nThresh = (int)((thresh*nCtg)/100);
					rs = s.executeQuery("select abs(" + colName +  ") from contig where pid > 0 order by " + colName + " asc limit " + nThresh);
					rs.last();
					double thresh2 = rs.getDouble(1);
					System.err.println(thresh + "% threshold:" + thresh2);
					thresh = thresh2;
				}
	        		
		        	int[] DE = new int[nCtg]; 
		        	rs = s.executeQuery("select abs(" + colName + ") from contig where pid > 0 order by ctgid asc");
		        	n = 0;
		        	while (rs.next())
		        	{
		        		double p = rs.getDouble(1);
		        		int de = (p <= thresh ? 1 : 0);
		        		DE[n] = de;
		        		n++;
		        	}
		        	re.assign("ctgDE",DE);
		        	
				cmd = "names(ctgDE) <- names";
				System.err.println("     " + cmd);
				re.eval(cmd);

				cmd = "gos <- vector(\"list\",nCtg)";
				System.err.println("     " + cmd);
				re.eval(cmd);

				cmd = "gonames <- names";
				System.err.println("     " + cmd);
				re.eval(cmd);

				String curID = "";
				Vector<Integer> ctgGOs = new Vector<Integer>();
				rs = s.executeQuery("select c.contigid,p.gonum from contig as c left join pja_unitrans_go as p " +
						" on p.ctgid=c.ctgid where c.pid > 0 order by c.ctgid asc");
				n = 1; // list index starts at one - why??
				while (rs.next())
				{
					String contigid = rs.getString(1);
					int gonum = rs.getInt(2);
					if (curID.equals("")) curID = contigid;
					if (!contigid.equals(curID))
					{
						if (ctgGOs.size() > 0)
						{
							StringBuilder str = new StringBuilder("");
							str.append("\"" + ctgGOs.firstElement() + "\"");
							for (int i = 1; i < ctgGOs.size(); i++)
							{
								str.append(",\"" + ctgGOs.get(i)  + "\"");
							}
							cmd = "gos[[" + n + "]] <- c(" + str.toString() + ")";
							//System.err.println("     " + cmd);
							re.eval(cmd);
							ctgGOs.clear();
						}
						else
						{
							cmd = "gonames[" + n + "] <- NA";
							re.eval(cmd);
							
						}
						curID = contigid;
						n++;
					}
					if (gonum > 0)	ctgGOs.add(rs.getInt(2));
				}

				cmd = "names(gos) = gonames"; 
				System.err.println("     " + cmd);
				re.eval(cmd);

				cmd = "np <- nullp(ctgDE,'','',ctgLens,FALSE)"; 
				System.err.println("     " + cmd);
				re.eval(cmd);
				
				cmd = "pvals <- goseq(np,'','',gos)";
				System.err.println("     " + cmd);
		        	RList results = re.eval(cmd).asList();
		        	String[] gonums = results.at("category").asStringArray();
		        	double[] pvals = results.at("over_represented_pvalue").asDoubleArray();
		        	assert(gonums.length == pvals.length);
		        	rs = s.executeQuery("show columns from go_info where field='" + colName + "'");
		        	if (!rs.first())
		        	{
		        		s.executeUpdate("alter table go_info add " + colName + " double default 0");	            
		        	}
		        	s.executeUpdate("update assem_msg set pja_msg=NULL");
		        	n = gonums.length;
		        	System.err.print("Uploading " + n + "     \r");
				PreparedStatement ps = conn.prepareStatement("update go_info set " + colName + "=? where gonum=?");
				for (int i = 0; i < gonums.length; i++)
				{
					int gonum = Integer.parseInt(gonums[i]);
					double score = pvals[i];
					ps.setDouble(1,score);
					ps.setInt(2,gonum);
					ps.execute();
					n--;
					if (n % 1000 == 0)
					{
						System.err.print("Uploading " + n + "     \r");
					}
				}
	        	}   
            System.err.println("Finished GOseq execution at: " + TimeHelpers.getElapsedTimeStr(startTime));
            System.err.println("\nThe console is in R, you may run R commands -- quit() when done, or perform another Execute.");
		}
		catch(Exception e)
		{
			System.err.println("Error - quit() out of R in terminal window.");
			ErrorReport.reportError(e, 
				"running GOseq -- may not have enough data to do anything meaningful.", true);
			return;
		}
	}
	
	/***********************************************
	 * JRI
	 */
	void initJRI()
	{
		try
		{
			if (re == null)
			{
				System.err.println("Startup R");
				Map<String, String> env = System.getenv();
				boolean flag=false;
			    for (String envName : env.keySet()) {
			            if (envName.contains("R_HOME")) {
			            		System.err.format("%s=%s\n\n", envName,env.get(envName));
			            		flag=true;
			            }
			    }
			    if (!flag) System.err.println("Warning: $R_HOME not set\n");
			    
				if (!Rengine.versionCheck()) {
				    System.err.println("** Version mismatch - Java files don't match library version.");
				    System.exit(1);
				}
				re= Rengine.getMainEngine();
				String[] args = {"--vanilla"};
				re = new Rengine(args, true, new TextConsole());
		        System.err.println("Rengine created, waiting for R");
		        
				// the engine creates R is a new thread, so we should wait until it's ready
		        if (!re.waitForR()) {
		            System.err.println("Cannot load R");
		            return;
		        }
			}
			else
			{
				System.err.println("Use existing R session");
				System.err.println("rm(list=ls())");
				re.eval("rm(list=ls())");
			}
		}
		catch(Exception e)
		{
			ErrorReport.die(e, "starting JRI");
		}
	}
	class TextConsole implements RMainLoopCallbacks
	{
	    public void rWriteConsole(Rengine re, String text, int oType) {
	        System.err.print(text);
	    }
	    
	    public void rBusy(Rengine re, int which) {
	        //System.err.println("rBusy("+which+")");
	    }
	    
	    public String rReadConsole(Rengine re, String prompt, int addToHistory) {
	        System.err.print(prompt);
	        try {
	            BufferedReader br=new BufferedReader(new InputStreamReader(System.in));
	            String s=br.readLine();
	            return (s==null||s.length()==0)?s:s+"\n";
	        } catch (Exception e) {
	            System.err.println("jriReadConsole exception: "+e.getMessage());
	        }
	        return null;
	    }
	    
	    public void rShowMessage(Rengine re, String message) {
	        System.err.println("rShowMessage \""+message+"\"");
	    }
		
	    public String rChooseFile(Rengine re, int newFile) {
			FileDialog fd = new FileDialog(new Frame(), (newFile==0)?"Select a file":"Select a new file", (newFile==0)?FileDialog.LOAD:FileDialog.SAVE);
			fd.setVisible(true);
			String res=null;
			if (fd.getDirectory()!=null) res=fd.getDirectory();
			if (fd.getFile()!=null) res=(res==null)?fd.getFile():(res+fd.getFile());
			return res;
	    }
	    
	    public void   rFlushConsole (Rengine re) {}
	    public void   rLoadHistory  (Rengine re, String filename) {}			
	    public void   rSaveHistory  (Rengine re, String filename) {}			
	}
	
	/** private variables **/
	private Connection conn;
	private SessionData selectedDB = null; 
	private static Rengine re = null;
	private HashSet<String> packages;
}
